The numerical simulation of physical phenomena is well known, and numerical solutions obtained from computer codes are used in many domains of industry. In the following example, we will recapitulate the steps which enter into the daily routine of the application of numerical simulations. The domain of application of the example is Computational Fluid Dynamics (CFD), but the methodology applies equally to other domains and applications. Consider a company which makes aircraft, and who wishes to design or to optimize a wing for the lift to drag ratio. The design department decides to use the tool of numerical simulations since wind tunnel tests are expensive. Their CFD department responds to this request in the following way.
1. A decision is made on which part of physics is important to describe in enough detail the phenomena in a meaningful simulation. In many fluid dynamics cases, the basic fluid dynamics equations like the Navier-Stokes equations are adequate. Perhaps a more simple description can be used. If viscous effects are less important, the Euler equations can be used. At low Mach number, an incompressible approximation can be used. Maybe the potential flow approximation is sufficient. On the other hand, a turbulence model may be needed, together with acoustic models, and a description of the deflection of the wing under aerodynamic forces.
The physical description is expressed in a mathematical model comprising of one or more equations which describe the relations between the physically relevant parameters (mass, density, energy, forces, impulse, pressure, volume, temperature, species concentration, charge density, magnetic field, . . . ). In general, the mathematical model uses spatial derivatives of one or more of the parameters, as well as time derivatives. Together with the number of space dimensions, the physical model determines the number of variables which are needed in the solution. In the case of the Navier-Stokes equations in three dimensions, five variables which depend on space and time describe the physical solution. These five variables can be e.g. the density, the momentum in the three directions and the energy. This set of variables is commonly known as conservative variables. Other combinations are possible, such as the density, the three velocity components and the pressure, which is known as the set of primitive variables. Yet another set of variables is known as entropy variables, since one of the variables is chosen to be proportional to the entropy. A transformation is always possible between different sets of independent variables.
2. The next step is to choose a computer program which can compute a solution to the problem which is compatible with the physical description of the problem. The simulation code will be run on a computer. While the physics and the mathematical model describe a continuum, the computer can only work with discrete values in a limited number of points. The solution of the mathematical model with a computer program therefore involves the procedure of discretization, which consists of two steps.
First, the continuous variable is replaced by a discrete variable in space. A grid defines points in space where the values of the parameters are stored.
The second step is to rewrite the relations between the continuous variables of the mathematical model in relations between discrete variables. In particular, the derivatives which appear in the mathematical model are replaced by their discrete equivalent, which are differences which are locally calculated using the values of the discrete variables of the surrounding grid points. This step involves an approximation and is not exact.
In general, a big enough company has the means and experience to develop and maintain an in-house computer code. Such a code is tuned to the specific needs of the company. It may compute solutions to only a very restricted number of physical problems, but it is likely to be up to date. The developments in the domain are followed, and e.g. new space discretizations, faster solution techniques such as accelerators and so on are coded. After testing, the code is ready for production runs. Smaller enterprises often can afford only a small CFD group, which uses programs from a specialized provider. This is normally a more general code, since it will have to accommodate a larger spectrum of clients with different applications.
3. The computation of the numerical solution is preceded by a pre-processing step. A geometrical description of the wing is used to generate a computational grid. This grid defines the points in space where the solution is calculated. Together with this grid an initial solution is generated. In the case of the Navier-Stokes equations, the five independent variables are given a certain value for each of the grid points, corresponding to e.g. uniform flow at the design Mach number. The grid and the initial solution are written to one or more files.
4. The next step is to run the simulation code. The program performs the following actions:                a) The program reads the grid, and stores for each point the coordinates in the memory of the computer.        b) It reads the initial solution and stores for each grid point the value of the physical parameters in memory.        c) Additional files are read which specify computational parameters such as the number of iterations or a convergence criterion, which physical model to use, which of the coded space discretizations to use, and additional information if necessary. Boundary conditions are specified which are used at the extremities of the part of space under consideration. They represent inflow and outflow from the domain, the way the wing influences the flow, and any other interaction of the flow with its surroundings.        d) Using the grid, the initial solution and the boundary conditions, it solves iteratively and approximately the discrete mathematical equations which describe the physical problem.        e) The code stops after a specified number of iterations, or when a solution is obtained which is accurate enough.        f) The code writes the solution to a file. This output solution consists of the value of physical variables in each of the grid points, where the number of physical variables corresponds to the problem. Another output file may be written which contains the grid.        
5. A post-processing program is used to read the grid file and the output file of step 4f. Normally this very general program permits to extract any desired physical parameter which can be computed from the solution. If the solution is written in the form of conservative variables, this program can compute the pressure, which is needed for the lift and drag. It can also compute variables like the Mach number, the temperature and so on.
6. The results are used to modify or improve the geometry.
Steps 3-6 are repeated until the design is optimal. In the case of insufficiency of the physical model, a more elaborate one is chosen, and coded if necessary. If there is insufficient resolution of the simulated physical phenomena on parts of the grid, a finer grid is generated.
The above describes the typical computation of a numerical solution and its use in an industrial environment. A similar procedure is followed when the purpose is to design or to optimize part of a car, a ship, a vehicle in general, a turbo machine, a ventilation system, a wind mill, a pump, a combustion engine, a channel, the mixture of steam and oil or water and oil in porous media for oil recovery, and so on, or to describe the evolution of flow such as in weather prediction, climate modeling or in a forest fire. In step 1, the physical model depends on the problem under consideration. In a magneto-hydrodynamics application, the description of the electromagnetic field enters. This entails the use of additional variables in the solution. In an application with different chemical species, the concentration of each of the species enters, together with a description of the transformation between the different chemical compounds. In multi-phase flow, such as e.g. oil and steam, or water and air, all fluids require equations and variables for a proper description.
The pre-processing and post-processing, steps 3 and 5 respectively, are normally separate programs. Occasionally, they can be incorporated in the numerical solver.
In step 3 the grid is generated once for a given geometry. In current practice a grid can contain many millions of points. Very large computational grids are normally composed of different blocks. Such a multi-block grid can then be distributed to multiple processors or even to different computers for a parallel computation. Some numerical simulation programs are able to refine or derefine automatically the grid locally using an error estimation based on the computed solution. Some programs use such a capability to compute a time-accurate simulation where a concentration of grid points follows a physical phenomena like e.g. a vortex or a shock. In the case of a time-accurate simulation, step 4f is repeated for each intermediate solution of interest. This is e.g. the case for the instationary behavior of a turbo machine where the radiated sound is computed. Another example is combustion.
The practical use of numerical simulations is illustrated in step 5, where the output of the program is transformed into a useful result. The output of the numerical simulation can serve many purposes, depending on the application. Quite often, this involves the pressure, since this relates to forces on the object under design. For a wing, this can be the Mach number for restricting the shock strength. For turbo machines, this can be the pressure on the blades, or the thermal load. For a heart valve, this can be the velocity field which needs to be smooth to avoid cluttering. For a forest fire simulation, this is again the velocity field to predict the likely progress of the fire.
Normally, the person who runs the program makes use of the output as illustrated in step 6. The output can also be fed back without human interference into the program. This is e.g. the case for inverse methods, where the desired pressure distribution is prescribed in the form of a goal. The geometry of the object and the related grid possess degrees of freedom which allow optimization to reach the goal. Remark that the purpose in the example above is the optimal geometry of the wing. The numerical solution is a convenient tool to achieve this goal.
At the core of the numerical simulation is the approximation of space derivatives as used in step 4d. The discrete mathematical equations which are used in that step are the subject of the invention.
For the computation of each value of a space derivative at a grid point, also called node, one of the methods, known under the name of Finite Differences, consists of the computation of an approximate value of each space derivative as a linear combination of values of the parameter(s) taken at each of the points of a set of grid points, called stencil, which are in general points surrounding the point of computation, and with weighting coefficients chosen in such a way that the best possible approximation results.
Another method is the Finite Volume method, which is very closely related to Finite Differences, and is based on an integral formulation using fluxes through a control volume. These fluxes are in turn a function of a number of the nodal unknowns, defining again a stencil. The number and relative importance of nodes involved in the fluxes allows for optimizing the approximation according to predefined criteria.
Yet another method is the Finite Element method, where the space is subdivided into elements, containing nodes at which the approximation of the physical parameter(s) u is stored. Basis functions 0, also called interpolation functions, are defined on the elements, which are used to approximate the physical value u on the element. The integral of the derivative is computed on the element, using a test function V, also called weighting function. The combination of X and V is then chosen in such a way that the approximation of the derivative is optimal.
Then, there are Distribution Methods. They share with Finite Elements the representation of the unknown(s) is on the element. One variant, the Residual Distribution method shares with the Finite Volume method a volume over which the integral of the derivative is computed. The integral of the derivative is split into parts which are distributed over the nodes. Another variant, the Flux Vector Distribution method uses like the Finite Volume method fluxes through a control volume. These fluxes are distributed over the nodes. In both cases, the characteristics and the optimization of the discretization follow from the choice of the distribution coefficients.
All these methods have in common that the approximation is ultimately dependent on the values of the parameter(s) at the nodes of the stencil, and that certain discretization parameters, which are inherent to the particular numerical framework under consideration, influence the accuracy of the approximation and the numerical properties of the numerical simulation of the physical phenomena.
The choice of the stencil, the values of the discretization parameters and the associated numerical method are commonly called the numerical scheme of the numerical simulation of physical phenomena.
Stencils are commonly classified according to the relative position of the points with respect to the advection direction in velocity fields. In an upwind discretization, the stencil has nodes which are upstream of the point where the derivative is computed. Another type of discretization is central, where the stencil is symmetric with respect to the advection speed for the point where the derivative is computed. Finally, there are many discretizations which are a combination of the two.
Concerning the grid, one can use a regular or structured mesh, formed by N families of rectilinear or curved lines, where the mesh points are located at the intersections. The number of families of lines corresponds to the number of space dimensions. Another type of grid, called unstructured grid, uses a connection table to identify the points of a stencil, and the points are in general less regularly positioned in space. However, unstructured grids can be locally perfectly regular, allowing to discern locally families of lines.
The accuracy of the simulation depends directly on the grid spacing. If the error in each space derivative reduces linearly with the mesh size, the scheme is said to be of first order. If the error reduces quadratically with the mesh size, the scheme is said to be of second order. In general, when the error reduces with the power M of the mesh size, the scheme is said to be of order M.
For a given application, the choice of the scheme is linked with practical constraints of the computation, called numerical constraints, which correspond to certain priorities imposed by the user on resulting errors of the simulation. Among the numerical constraints, one may refer to consistency, convergence and stability of the simulation, the desired order of the scheme (industrial applications require an order of at least 2), the diffusion and dispersion properties of the numerical scheme, and computational cost (which may impose a maximum size of the stencil, generally indicated by compactness).
The order of the approximation and the extent of the stencil are related. On a stencil of given extent, the order of accuracy of an approximation is limited. In other words, for an approximation of a given order, the stencil has to have a certain minimum extent. It depends on the numerical constraints of the application if the size of the stencil or the order of the approximation is predominant. When the appropriate combination of the size of the stencil and the desired accuracy is chosen, remaining degrees of freedom can be used for optimizing other numerical constraints.
An extended description of the above, and of numerical methods in general can be found in (D1), “Numerical computation of internal and external flow”, C. Hirsch, J. Wiley & Sons, New York, Vol. 1, Fundamentals of numerical discretization, 1988 and Vol. 2, Computational methods for inviscid and viscous flows, 1990.
There are various methods to obtain the discretization parameters appearing in the specific numerical framework. In general, a truncated Taylor series expansion with respect to the computational point is used for the value of the physical parameter(s). Given that only the first terms of the Taylor series are used, the value obtained is accompanied with an error, called truncation error, which depends in particular on the degree of the last retained term in the series. The discretization parameters are then chosen such that they optimize the error in the discretization of the derivative, which is a function of the values resulting from the truncated Taylor series. Other criteria follow from considering the Fourier transform of the discretization, thereby choosing the discretization parameters in such a way that the behavior of certain Fourier components is optimized. This means tuning the dissipation and dispersion of the numerical scheme, optimizing for the stability of the numerical scheme, or for the damping of certain components resulting in accelerated convergence. Other optimizations related to the choice of the discretization involve computational cost, ease of use in multi-block codes, ease of coding and others.
As pointed out, for each application, the most appropriate choice has to be made on:                the numerical formulation with associated discretization parameters which is used for the simulation,        the stencil itself, and notably the number of points and their relative position with respect to the point of computation,        the optimization of the discretization, e.g. based on the truncation error in the truncated Taylor series expansion used in the computation of the discretization parameters, or based on the behavior of certain Fourier components, or other.        
Many methods for numerical simulation have been proposed and are used in industry with more or less satisfactory performances. The above mentioned books D1 give an overview of the mathematical models, discretization techniques, numerical schemes and solution methods for the numerical simulation of space derivatives in flow phenomena. In particular, in D1 Volume 1, part 2, chapter IV, pages 167-180, the fundamentals are described of the method of Finite Differences, on a regular grid in one dimension. As explained, the weighting coefficients used in the linear combination can be determined if the number of points of the stencil is coherent with the order of the desired discretization. In particular, the book mentions the general method of Hildebrand for determining Finite Difference formulas.
An example is a general one-dimensional discretization for a derivative at node i which uses the stencil between the nodes i−m and i+n, where m and n are given positive integers. On the stencil S=(i−m,i−m+1, . . . , i−1,i,i+1, . . . , i+n−1,i+n), the approximation of the first derivative can be written as
                                                                        u                x                            =                            ⁢                                                1                                      Δ                    ⁢                                                                                  ⁢                    x                                                  ⁢                                  (                                                                                    a                                                  -                          m                                                                    ⁢                                              u                                                  i                          -                          m                                                                                      +                                                                  a                                                                              -                            m                                                    +                          1                                                                    ⁢                                              u                                                  i                          -                          m                          +                          1                                                                                      +                    …                    +                                                                                                                                        ⁢                                                                    a                                          -                      1                                                        ⁢                                      u                                          i                      -                      1                                                                      +                                                      a                    0                                    ⁢                                      u                    i                                                  +                                                      a                    1                                    ⁢                                      u                                          i                      +                      1                                                                      +                …                +                                                                                                                        ⁢                                                                            a                                              n                        -                        1                                                              ⁢                                          u                                              i                        +                        n                        -                        1                                                                              +                                                            a                      n                                        ⁢                                          u                                              i                        +                        n                                                                                            )                            .                                                          (        1        )            
The coefficients aj, j∈S determine the numerical properties of the discretization. A general description of the above discretization can be found in D1 which discusses the method for obtaining the coefficients aj. Expression (1) also generalizes to higher derivatives.
The conjunction in two or more dimensions of the above one-dimensional discretization has a small error if the preferential direction points from one grid point to another, following a grid line. In between, the error becomes quite large. The reason is that a derivative which is related to a preferential direction has to be approximated with grid-based derivatives. Grid-based derivatives are using values of the physical unknown(s) which are stored at predefined positions of the grid. A conjunction of one-dimensional discretizations uses only points on the grid lines, and not in between.
This is illustrated in FIG. 1, where a part of a two-dimensional structured grid is shown. The indices are for the x-direction 0,1, . . . , i−1, i, i+1, . . . , imax and for the y-direction 0,1, . . . , j−1, j, j+1, . . . , jmax, and the nodal unknowns are indicated by ui,j for node i, j. Indicated are two directions, {right arrow over (a)}1 and {right arrow over (a)}2. Both point from one grid point to another, but only {right arrow over (a)}1 follows a grid line, and the discretization has a small error. For the direction of {right arrow over (a)}2, the error is maximal when using a conjunction of one-dimensional discretizations.
The report (D2), “Progress in multidimensional upwind differencing”, B. van Leer, NASA Contractor Report 189708, ICASE report #92-43, September 1992, NASA Langley Research Center, Hampton, Va., describes different known approaches for the simulation of convection phenomena in two dimensions which go beyond one-dimensional methods. In this framework, different state of the art approaches have been proposed to optimize the methods of numerical simulation in two dimensions, with more or less satisfactory performances. Only the methods using an unstructured grid have been exploited in practice using a low order upwind scheme.
In the majority of practical applications, the physical phenomena take place in three dimensional space, which makes it indispensable to realize a numerical simulation of at least three dimensions. In practice, it may be necessary to use more than three dimensions, e.g. to incorporate time as additional dimension, taking into account the fact that an unsteady physical process in N dimensions can be modeled by a N+1 dimensional steady process.
However, actually all the industrial methods of numerical simulation in three dimensions (cf. e.g. the documents D1 and D2) are limited to a conjunction of one-dimensional methods.
Besides these one-dimensional methods, certain attempts of truly multi-dimensional numerical simulations in three dimensions have been reported. Nevertheless, taking into account the complexity of the problem presented, and the difficulties encountered in two dimensions, the use of a truly multi-dimensional method in three or more dimensions is considered to be extremely complex.
The paper (D3) “Optimum positive linear schemes for advection in two and three dimensions”, P. L. Roe and D. Sidilkover, SIAM J. Num. Anal., December 1992, #6, vol. 29, pages 1542-1568, constitutes one of the rare publications which treat until now a scheme for the simulation of advection phenomena, considering a class of first order upwind schemes. In the case of three dimensions, the authors of the paper, like all authors until now, consider it suitable to choose on a structured grid a minimal upwind stencil formed by a cube of eight points.
Having chosen and predefined this stencil, the different weighting coefficients are computed to obtain the best simulation possible, which is the case of the scheme called “N scheme”. Denote the components of the vector {right arrow over (a)} of the preferential direction in N dimensions by a, b, c, . . . , that is {right arrow over (a)}=(a, b, c, . . . )T. One-dimensional upwind schemes take the discretization stencil based on a≧0 or a≦0. Applying a conjunction of one-dimensional discretization in N dimensions means choosing the stencil depending on the sign of a, b, c, . . . separately. The N scheme uses instead the combination a≧b≧c≧ . . . ≧0 and points of the stencil which do not lie on a grid line through the point of computation P. Nevertheless, the results obtained with the scheme given in the paper are largely insufficient in practice. Indeed, one cannot obtain a high order scheme with this method on such a small stencil. Furthermore, the simulation in three dimensions generates specific problems which are not encountered in two dimensions.
Moreover, since the date of this paper, no author has effectively imagined that it is possible to obtain higher order accurate discretizations using a larger stencil which includes points beyond the grid lines on a structured grid in three dimensions.
However, as indicated by the publication (D4) “Numerical Aerodynamics: Past Successes and Future Challenges from an Industrial Point of View”, David Paul Hills, Computational Methods in Applied Sciences, 1996, ECCOMAS, pages 166-173, one-dimensional methods as used in industry have reached their limits, and real gains should be sought in multi-dimensional methods. In spite of continued research since 1992, no such method has been found yet.
Furthermore, the majority of the publications on this subject since 1990 indicate, like the previous, that the true multi-dimensional methods imply upwind discretizations on unstructured grids. Neither multi-dimensional efforts on structured grids nor attempts to find central space discretizations with a reduced error have resulted in numerical simulations which can be applied in an industrial environment.
The invention aims to dispose of these shortcomings by proposing a simulation method which provides a simulated numerical representation in N dimensions, with N≧3, of spatial pth derivatives Dp, p≧1, of at least one physical parameter u of a phenomenon to which is associated a velocity field d which determines a preferential direction in every point of at least a part of space, and which permits                to reduce significantly the error in the numerical simulation,        a consistent, stable and convergent simulation,        a locally monotone simulation,        the simulation of spatial pth derivatives, with p>1, with an error of a desired order, incorporating numerical constraints as desired.        
The invention aims also to propose a method which                is compatible with current technologies, notably computation time and memory requirements,        is compatible with the simulation of nonlinear phenomena,        can be used for systems of equations,        can be used with acceleration techniques like e.g. multi-grid, GMRES and preconditioning,        and which allows simulations in Finite Element formulations, Finite Volume formulations, Finite Difference formulations, Residual Distribution formulations, Flux Vector Distribution formulation, or others.        
The invention also aims at providing a global simulation method of a physical phenomenon and a data processing system, a software, in particular in the form of a digital storage medium, to implement such methods.